# Metview Macro

# **************************** LICENSE START ***********************************
#
# Copyright 2012 ECMWF. This software is distributed under the terms
# of the Apache License version 2.0. In applying this license, ECMWF does not
# waive the privileges and immunities granted to it by virtue of its status as
# an Intergovernmental Organization or submit itself to any jurisdiction.
#
# ***************************** LICENSE END ************************************


# This macro computes and plots the difference between a forecast
# and an analysis field, specified by the user via a user interface
# Run the macro and the interface will pop up.


# define the user interface components
a = any(
        name    : "base date",
        default : "-2"
        )

b = option_menu(
        name    : "base time",
        values  : ["0", "12"],
        default : "12"
        )

c = option_menu(
        name    : "step (hours)",
        values  : ["0", "12", "24", "36"],
        default : "24"
        )

d = option_menu(
        name        : "parameter",
        values      : ["UV", "t", "z"],
        default     : "t"
        )

e = any(
        name        : "area",
        help        : "help_input",
        input_type  : "area",
        default     : [30,-25,50,65]  # S, W, N, E
        )


# retrieve the input values
input = dialog([a, b, c, d, e])

if input <> nil then
    base_date   = input["base date"]
    base_time   = input["base time"]
    step        = input["step (hours)"]
    par         = input["parameter"]
    the_area    = input["area"]

    # translate "UV" into a list, since that is what 'retrieve' understands
    if (par = "UV") then
        par = ["u", "v"]
    end if

else
    fail("macro failed to get input elements")
end if


# Define some visual definitions

pos = mcont(
        contour_line_thickness       : 2,
        contour_line_colour          : "red",
        contour_highlight            : "off",
        contour_level_selection_type : "level_list",
        contour_max_level            : 10,
        contour_min_level            : 0.5,
        contour_level_list           : [0.5,1,2,4,10]
        )

neg = mcont(
        contour_line_thickness       : 2,
        contour_highlight            : "off",
        contour_level_selection_type : "level_list",
        contour_max_level            : -0.5,
        contour_min_level            : -10,
        contour_level_list           : [-10,-4,-2,-1,-0.5]
        )


coloured_wind = mwind(
        legend                                : "on",
        wind_arrow_legend_text                : "m/s",
        wind_advanced_method                  : "on",
        wind_advanced_colour_max_level_colour : "red",
        wind_advanced_colour_min_level_colour : "blue",
        wind_advanced_colour_direction        : "clockwise",
        wind_thinning_factor                  : 1
        )

acoast = mcoast(
        map_coastline_resolution        : "low",
        map_coastline_thickness         : 3,
        map_coastline_land_shade        : "on",
        map_coastline_land_shade_colour : "grey",
        map_coastline_sea_shade         : "on",
        map_coastline_sea_shade_colour  : "RGB(0.9,0.95,1)",
        map_grid_longitude_increment    : 10
        )

ps_atlantic = geoview(
        map_projection      : "polar_stereographic",
        map_area_definition : "corners",
        area                : the_area,
        coastlines          : acoast
        )


# Define the display window

page = plot_page(
        view  : ps_atlantic
        )

dw = plot_superpage(
        custom_width    : 29.7,
        custom_height   : 21.0,
        pages           : page
        )


# Check which visual definition to use: contour or wind arrows?

if (par = ['u', 'v']) then
    visdef = coloured_wind
else
    visdef = [pos, neg]
end if

# retrieve the data

fc_date = date(base_date) + hour(base_time) + hour(step)
print(fc_date)

an_data = retrieve
(
    date: yyyymmdd(fc_date),
    time: hour(fc_date)*100,
    parameter: par,
    grid: [1,1]
)

fc_data = retrieve
(
    type: 'fc',
    step: step,
    date: base_date,
    time: base_time*100,
    parameter: par,
    grid: [1,1]
)



# Derive the difference field

fa_diff = fc_data - an_data


# Plot the result

plot (dw, fa_diff, visdef)


